Bare Necessities

Author

Adam Kemberling

Published

November 27, 2024

Code
# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
library(nlme)
library(ggimage)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "white", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom"))



# labels for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Flag effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * area * season) ) %>% 
    mutate(
      area = factor(group, levels = area_levels_long_all),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

Everything in this markdown uses the size spectra exponent values before I shifted the wmin parameter. The wmin values are fixed at 1g across all groups.

Code
# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>% left_join(area_df)

# Load the three regional datasets
wigley_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_bodymass_spectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))

# Load Shelfwide
wigley_bmspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))



# Join them together
wigley_bmspectra_all <- bind_rows(wigley_bmspectra, wigley_bmspectra_shelf) %>% 
  mutate(metric = "bodymass_spectra",
         community = "Wigley Subset",
         survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  left_join(area_df) %>% 
  mutate(area = factor(area, area_levels_long_all))


# both metrics+communities
spectra_results <-  wigley_bmspectra_all %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels_all),
    metric = fct_rev(metric))
Code
# # Model Datasets for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# load shelfwide stuff too
wigley_bmspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))

# Bottom Temperatures
bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))


# Bolt that on
wigley_bmspectra_model_df <- bind_rows(
  wigley_bmspectra_model_df,
  left_join(
    wigley_bmspectra_shelf, 
    bot_temps, 
    by  = join_by("est_year" == "year", survey_area, season))) %>% 
  left_join(area_df) %>% 
  mutate(area = if_else(survey_area == "Northeast Shelf", "Northeast Shelf", area),
         area = factor(area, levels = area_levels_long_all),
  season = fct_rev(season))
Code
# From Andy Beet - Crabs removed
landings_annual <- read_csv(here::here("Data/processed/BEET_GARFO_regional_finfish_landings.csv")) %>% 
  select(survey_area,
         est_year = year,
         total_weight_lb = total_live_lb)

# Get overall landings
shelf_landings_annual <-  mutate(landings_annual, survey_area = "Northeast Shelf") %>% 
    group_by(est_year, survey_area) %>% 
    summarise(total_weight_lb = sum(total_weight_lb, na.rm = T)) %>% 
    mutate(area = survey_area)

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

Temperature and Size Structure

There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.

  • Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
  • James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes

The metabolic theory of ecology provides a basis for “why” the above rules might be observed.

  • Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.

Seasonal Movements and Metabolic Efficiency: MTE details how metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.

This theory would suggest that larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and to reduce the energetic costs of thermoregulation in warmer periods.

From what I see at first-glance the size distribution along the shelf is responding consistently with temperature changes (Spring to Fall) on shorter time scales than hypotheses that predict temperature related changes in growth.

Code
wigley_bmspectra_model_df %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b, color = season)) +
  geom_point(aes(shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  ggforce::geom_mark_ellipse(aes(fill = season), alpha = 0.1) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  guides(
    shape = guide_legend(nrow = 2),
    fill = guide_legend(nrow = 2),
    color = guide_legend(nrow = 2)) +
  theme(legend.box = "horizontal") +
  labs(
    title = "Bergmann's Rule + MTE",
    subtitle = "The broad relationship between temperature and size distribution is bolstered by\nseasonal community changes in alignment with seasonal temperature differences.",
    shape = "Survey Area",
    fill = "Season",
    color = "Season",
    x = "Bottom Temperature")

Importance of Seasonality

The extent that the regional spectra trends follow a negative linear relationship to increasing temperature is largely due to seasonal differences.

We know that individuals may migrate and/or occupy different habitats seasonally, leading to situations where the community compositions from spring to fall are an apples to oranges comparison. Seasonal movements may be be partially explained by expectations from the MTE, but they also could be separate ecological processes like seasonal foraging or spawning opportunities.

Either way, seasonal differences are large and they justify treating survey seasons as separate groups with independent trends.

Anomalies Model

I wanted to try and avoid conflating that a relationship between spectra and bottom temperature means that “warming” has altered the size spectra. The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.

This is the model I would suggest gets the most at “Warming” impacting regional size distributions:

Fixed effects:

  • Season
  • Survey area
  • Bottom temperature anomaly

Error structure: AR1 autocorrelative errors

By using the season * region interaction groups, we can characterize the typical bottom temperature for each combinations using the long-term seasonal average, and quantify warming using departures from those long-term averages.

The model looks like this as an equation (90% sure), where \(\alpha\) is the intercept, \(\beta x_{tjs}\) is our coefficient for a continuous independent variable (ex. bottom temperature), and \(e_{tjs\)` is our error:

\[y_{tjs} = \alpha + \beta x_{tjs} + e_{tjs}\] and our error depends on the previous timestep:

\[e_{tjs} = \phi e_{t-1,js} + w_{tjs}\] With \(w_{tjs}\) adding noise.

\[w_{tjs} \sim N(0, \sigma^2) \] All of this happening for each \(year_t\) within \(region_j\) & \(season_s\).

I need to consult the Zuur paper about equations in ecology again, but this seems close.

Code
# Just drop the random effect part for now
anom_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp_anom), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(anom_model_ar)$corStruct


# # check the model
# check_model(anom_model_ar)
# plot(check_collinearity(anom_model_ar))
# plot(
#   check_collinearity(
#     nlme::gls(
#       b ~ area + season + scale(bot_temp_anom), #+ scale(log(total_weight_lb)), 
#       data = wigley_bmspectra_model_df, 
#       correlation = corAR1(form = ~ est_year | area/season)
# ))
#   )


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    anom_model_ar, 
    terms = ~ bot_temp_anom * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = anom_model_ar,
  specs =  ~ area * season,
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_anoms <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Based on this model I would say there is evidence of a warming effect in:
1. GOM-Spring
2. GB-Spring
3. GB-Fall
4. SNE-Fall

Bottom Temperature Model

We were originally using absolute temperatures, which are scaled, and I was curious with the AR error structure how different this would be than the anomaly model.

I just had to see if the same model is better performing or leads to different results when it knows absolute temperatures.

Code
# Just drop the random effect part for now
temp_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct
# 
# # check the model
# check_model(temp_model_ar)

# # Collinearity without interactions
# plot(
#   check_collinearity(
#      nlme::gls(
#         b ~ area + season + scale(bot_temp), 
#         data = wigley_bmspectra_model_df, 
#         correlation = corAR1(form = ~ est_year | area/season))
# ))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    temp_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = temp_model_ar,
  specs =  ~ area * season,
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_results <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Seasonal Bottom Temperature")
local_btemp_results

Verdict: Temperature itself performs better, and is simpler to produce/explain.

Code
data.frame(
  "Model" = c("Temperature Model", "Anomalies Model"),
  "AIC" = c(AIC(temp_model_ar), AIC(anom_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()
Model AIC Delta-AIC
Temperature Model -1206.481 0.00000
Anomalies Model -1184.586 21.89431

Thoughts:

Consistent seasonal shifts in the Shelf-wide spectra suggests to me that the community is more responsive to within-year time scale factors and maybe less sensistive to long-term trends in temperature than I originally expected. This also raises questions about how the spring/fall size distribution differences are being achieved and whether to view them as expressions of the same “community” or separately.

From the decadal variability work we have reason to believe that individual species are both not (perfectly) tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

In light of that knowledge, there must be some compensatory mechanisms happening across this broader community for the temperature+spectra relationship to be so limited/situational.

Region * Season Pairwise Posthoc

The following figures and tables compare marginal means and post-hoc comparisons for the model above using bottom temperature.

What are the differences across spatial scales and seasonally?

Code
# Region and seasonal posthocs
regseas_phoc <- emmeans(
  temp_model_ar, 
  list(pairwise ~ area + season), 
  adjust = "tukey",
  type = "response")



# Custom Plot
(intercepts_emmeans <- regseas_phoc$`emmeans of area, season` %>% 
    as_tibble() %>%
    ggplot() +
    geom_pointrange(
      aes(area, emmean, ymin = lower.CL, ymax = upper.CL, color = season),
      position = position_dodge(width = 0.25), alpha = 1) +
    #coord_flip() +
    scale_color_gmri() +
    labs(
      y = "Marginal Mean: Mass Spectra Slope",
      x = "Area",
      title = "Wigley Species",
      color = "Season",
      subtitle = "Group marginal means plot"))

Code
# Table for pairwise test results
regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  arrange(p.value) %>% 
  gt() %>% 
  tab_header(title = "Factor Groups Pairwise Post-Hoc") %>% 
  tab_style(
    style = list(
      cell_fill(color = "#9AFF9A"), 
      cell_text(color = "black")),
    locations = cells_body(
      rows = p.value <= 0.05) 
  )
Factor Groups Pairwise Post-Hoc
Contrast estimate SE df t.ratio p.value
Southern New England Spring - Northeast Shelf Fall 0.139 0.030 469.383 4.715 0.000
Southern New England Spring - (Mid-Atlantic Bight Fall) 0.237 0.045 476.537 5.263 0.000
(Mid-Atlantic Bight Spring) - Northeast Shelf Fall 0.157 0.023 455.375 6.921 0.000
(Mid-Atlantic Bight Spring) - (Mid-Atlantic Bight Fall) 0.255 0.041 476.029 6.240 0.000
Northeast Shelf Fall - Gulf of Maine Fall -0.143 0.029 465.859 -5.004 0.000
Northeast Shelf Fall - Georges Bank Fall -0.182 0.031 471.728 -5.895 0.000
Gulf of Maine Fall - (Mid-Atlantic Bight Fall) 0.240 0.044 444.388 5.418 0.000
Georges Bank Fall - (Mid-Atlantic Bight Fall) 0.279 0.046 477.818 6.089 0.000
Northeast Shelf Spring - (Mid-Atlantic Bight Fall) 0.201 0.049 473.748 4.104 0.002
Georges Bank Fall - Southern New England Fall 0.121 0.036 477.467 3.412 0.024
Gulf of Maine Spring - (Mid-Atlantic Bight Fall) 0.170 0.050 476.620 3.391 0.026
(Mid-Atlantic Bight Spring) - Southern New England Fall 0.097 0.029 475.519 3.367 0.028
Southern New England Fall - (Mid-Atlantic Bight Fall) 0.158 0.047 415.081 3.335 0.031
Northeast Shelf Spring - Northeast Shelf Fall 0.103 0.035 367.778 2.934 0.100
Georges Bank Spring - (Mid-Atlantic Bight Fall) 0.158 0.057 472.321 2.800 0.139
Gulf of Maine Spring - Georges Bank Fall -0.109 0.039 477.455 -2.757 0.154
Georges Bank Spring - Georges Bank Fall -0.121 0.047 476.957 -2.559 0.240
Gulf of Maine Spring - (Mid-Atlantic Bight Spring) -0.084 0.033 472.244 -2.523 0.259
Gulf of Maine Fall - Southern New England Fall 0.082 0.034 472.581 2.454 0.297
Southern New England Spring - Southern New England Fall 0.079 0.034 476.466 2.294 0.395
Georges Bank Spring - (Mid-Atlantic Bight Spring) -0.096 0.042 477.515 -2.280 0.404
Northeast Shelf Fall - (Mid-Atlantic Bight Fall) 0.097 0.044 474.070 2.219 0.445
Northeast Shelf Spring - Georges Bank Fall -0.078 0.038 471.924 -2.073 0.548
Gulf of Maine Spring - Northeast Shelf Fall 0.073 0.037 473.955 1.970 0.621
Gulf of Maine Spring - Gulf of Maine Fall -0.070 0.038 476.888 -1.853 0.701
Northeast Shelf Fall - Southern New England Fall -0.060 0.033 471.108 -1.832 0.715
Georges Bank Spring - Gulf of Maine Fall -0.082 0.046 477.356 -1.790 0.741
Gulf of Maine Spring - Southern New England Spring -0.066 0.038 477.822 -1.729 0.779
Northeast Shelf Spring - (Mid-Atlantic Bight Spring) -0.054 0.031 464.727 -1.716 0.786
Georges Bank Spring - Southern New England Spring -0.078 0.046 477.143 -1.694 0.799
Georges Bank Spring - Northeast Shelf Fall 0.061 0.045 474.606 1.348 0.942
Southern New England Spring - Georges Bank Fall -0.042 0.032 476.714 -1.301 0.953
Gulf of Maine Fall - Georges Bank Fall -0.039 0.031 440.843 -1.235 0.966
Northeast Shelf Spring - Gulf of Maine Fall -0.039 0.036 468.487 -1.095 0.985
Northeast Shelf Spring - Southern New England Fall 0.043 0.039 471.482 1.092 0.985
Northeast Shelf Spring - Southern New England Spring -0.036 0.037 470.472 -0.980 0.993
(Mid-Atlantic Bight Spring) - Georges Bank Fall -0.024 0.026 477.532 -0.918 0.996
Northeast Shelf Spring - Georges Bank Spring 0.043 0.050 474.220 0.848 0.998
Northeast Shelf Spring - Gulf of Maine Spring 0.030 0.043 473.560 0.709 0.999
Southern New England Spring - (Mid-Atlantic Bight Spring) -0.018 0.025 434.985 -0.723 0.999
Gulf of Maine Spring - Georges Bank Spring 0.012 0.051 422.378 0.235 1.000
Gulf of Maine Spring - Southern New England Fall 0.013 0.041 474.664 0.308 1.000
Georges Bank Spring - Southern New England Fall 0.001 0.049 437.668 0.012 1.000
Southern New England Spring - Gulf of Maine Fall -0.003 0.030 477.749 -0.109 1.000
(Mid-Atlantic Bight Spring) - Gulf of Maine Fall 0.015 0.024 472.981 0.622 1.000
Code
regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  filter(p.value < 0.05)  %>% 
  separate(Contrast, into = c("group1", "group2"), sep = " - ") %>% 
  mutate(
    group1 = str_remove_all(group1, "[(]|[)]"),
    group2 = str_remove_all(group2, "[(]|[)]"),
    statement = 
      if_else(estimate > 0, " is shallower than ", " is steeper than ")) %>% 
  select(group1, statement, group2) %>% 
  arrange(statement, group1) %>% 
  gt()
group1 statement group2
Georges Bank Fall is shallower than Southern New England Fall
Georges Bank Fall is shallower than Mid-Atlantic Bight Fall
Gulf of Maine Fall is shallower than Mid-Atlantic Bight Fall
Gulf of Maine Spring is shallower than Mid-Atlantic Bight Fall
Mid-Atlantic Bight Spring is shallower than Northeast Shelf Fall
Mid-Atlantic Bight Spring is shallower than Southern New England Fall
Mid-Atlantic Bight Spring is shallower than Mid-Atlantic Bight Fall
Northeast Shelf Spring is shallower than Mid-Atlantic Bight Fall
Southern New England Fall is shallower than Mid-Atlantic Bight Fall
Southern New England Spring is shallower than Northeast Shelf Fall
Southern New England Spring is shallower than Mid-Atlantic Bight Fall
Northeast Shelf Fall is steeper than Gulf of Maine Fall
Northeast Shelf Fall is steeper than Georges Bank Fall

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

Code
# Plot annual totals
bind_rows(list(
  shelf_landings_annual,
  landings_annual)) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_long_all)) %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.8) +
  geom_ma(
    aes(color = "5-Year Moving Average"), 
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries", linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Landings Autoregressive model

I wanted to follow on the bottom temperature anomalies model here for landings:

Code
# Just drop the random effect part for now
f_model_ar <- nlme::gls(
  b ~ area * season + area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
    dplyr::filter(est_year > 1978),
  correlation = corAR1(form = ~ est_year | area/season)
)



# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    f_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = f_model_ar,
  ~area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite")


# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = f_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(mod2_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free") +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Code
# temp model is better
# AIC(temp_model_ar)
# AIC(f_model_ar)

This is our old friend “higher landings when spectra was shallower” that we’ve seen before.

Temperature and Landings Together

If we want to just throw them both in because we’ve done what we can and want to hold them both in the frame at once we get this:

Code
# Just drop the random effect part for now
full_model_ar <- nlme::gls(
  # Model
  model = b ~ area * season * scale(bot_temp) + area * scale(log(total_weight_lb)), 
  
  # Data
  data = wigley_bmspectra_model_df %>%
    drop_na(total_weight_lb) %>% 
    dplyr::filter(area != "Northeast Shelf"),
  
  # Autocorrelation
  correlation = corAR1(form = ~ est_year | area/season))

# # Collinearity without interactions
# plot(
#   check_collinearity(
#     nlme::gls(
#         b ~ area + season + scale(bot_temp) +  scale(log_land),
#         data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%
#                  filter(area != "Northeast Shelf") %>% 
#           mutate(log_land = log(total_weight_lb)),
#         correlation = corAR1(form = ~ est_year | survey_area/season)
#       )
#     )
#   )


# # confidence intervals on phi
# intervals(full_model_ar)$corStruct


# temp model is much better
# AIC(temp_model_ar)
# AIC(full_model_ar)

Before proceeding it is worth noting that this model is less parsimonious than a temperature or anomaly only model:

Code
data.frame(
  "Model" = c(
    "Temperature Model", 
    "Anomalies Model", 
    "Landings Model", 
    "Landings and Temperature"), 
  "AIC" = c(
    AIC(temp_model_ar), 
    AIC(anom_model_ar),
    AIC(f_model_ar),
    AIC(full_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()
Model AIC Delta-AIC
Temperature Model -1206.4805 0.00000
Anomalies Model -1184.5862 21.89431
Landings Model -744.1868 462.29373
Landings and Temperature -897.3938 309.08670
Code
# Clean up the plot:
# Plot the predictions over data
temp_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
temp_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area * season,
  mode = "appx-satterthwaite",
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  filter(area != "Northeast Shelf") %>% 
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")



# Plot over observed data
local_btemp <- temp_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(temp_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area!= "Northeast Shelf"),
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp

Code
# Clean up the plot:
# Plot the predictions over data
f_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
f_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite"
  ) %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- f_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(f_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area != "Northeast Shelf"),
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

Code
zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()

Cold Pool Index

There are several cold pool indices available that might be helpful in adding context to changes in MAB. Below are the Dupontavice et al. 2022 methods that are used in the SOE report.

The Cold Pool Index (Model_CPI) was adapted from Miller et al. (2016). Residual temperature was calculated in each grid cell, \(i\), in the Cold Pool domain as the difference between the average bottom temperature at the year \(y\) (\(T_{i,y}\)) and the average bottom temperature over the period 1972–2019 (\(\bar{T}_{i, 1972-2019}\)⁠⁠) between June and September. \(CPI\) was calculated as the mean residual temperature over the Cold Pool domain such that

\[CPI_y~=~\frac{\sum_{i=1}^{n}~(T_{i,y}-\bar{T}_{i,~1972-2019})}{n}\]

where n is the number of grid cells over the Cold Pool domain.

Takeaway:

\(CPI\) was calculated as the mean residual temperature over the Cold Pool domain

From this index we can see that the temperatures over the cold pool domain (which include the MAB) have increased over time.

Code
cpool <- ecodata::cold_pool
cpool %>% 
  filter(!str_detect(Var, "se_")) %>% 
  ggplot(aes(Time, Value)) +
  geom_line() +
  facet_grid(Var~EPU, scales = "free")

Code
# # There's an SF
# cpool_sf <- ecodata::cold_pool_sf %>% 
#   mutate(id = row_number())
# 
# ggplot(cpool_sf) + geom_sf(aes(fill = id), alpha = 0.2)

Driver Correlations